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We determine the speed of a crystallisation (or more generally, a solidification) front as it advances 
into the uniform liquid phase after the system has been quenched into the crystalline region of the 
phase diagram. We calculate the front speed by assuming a dynamical density functional theory 
model for the system and applying a marginal stability criterion. Our results also apply to phase 
field crystal (PFC) models of solidification. As the solidification front advances into the unstable 
liquid phase, the density profile behind the advancing front develops density modulations and the 
wavelength of these modulations is a dynamically chosen quantity. For shallow quenches, the selected 
wavelength is precisely that of the crystalline phase and so well-ordered crystalline states are formed. 
However, when the system is deeply quenched, we find that this wavelength can be quite different 
from that of the crystal, so that the solidification front naturally generates disorder in the system. 
Significant rearrangement and ageing must subsequently occur for the system to form the regular 
well-ordered crystal that corresponds to the free energy minimum. Additional disorder is introduced 
whenever a front develops from random initial conditions. We illustrate these findings with results 
obtained from the PFC. 



I. INTRODUCTION 

It is important to understand the formation kinetics 
of a solid from the liquid phase when it is cooled be- 
low its freezing temperature Tf, because the microscopic 
structure of the solid can depend strongly on the forma- 
tion pathway. If the liquid is only slightly cooled below 
Tf, then the solid forms via nucleation and growth [I], 
generally leading to well-ordered crystalline solids. De- 
pending on the material and the degree of cooling below 
Tf, the formation of dendrites and other complex mi- 
crostructures is possible [2rllj • When the liquid is rapidly 
quenched (supercooled) to a temperature sufficiently far 
below Tf there is no nucleation barrier against the liq- 
uid forming a solid. In this situation, the solid that is 
formed can be amorphous, with little or no long-range 
order, rather than a regular ordered crystal. 

Classical density functional theory (DFT) |6] is a 
widely used microscopic theory capable of describing 
equilibrium aspects of melting, freezing and also the in- 
terfaces between the liquid and solid phases [5]. In 
conjunction with the recent development of a dynami- 
cal density functional theory (DDFT) [5HP2]. which is a 
theory that requires as input the free energy functionals 
from equilibrium DFT, these theories have been shown 
to be able to describe the dynamics of crystal forma- 
tion [T3]. A related approach, that has been developed 
and studied extensively over the last decade or so, is 
the phase field crystal (PFC) H IT3H2!?] approach for 
modeling the atomic structure of crystalline materials. 
The PFC may be derived from the DDFT by assum- 
ing a (local) gradient expansion approximation for the 
Helmholtz free energy functional for the system and lin- 
earizing the density-dependent mobility pre-factor in the 
DDFT equation [13 HI]. DFT and DDFT are theories 



for the one body density distribution p(x) of the particles 
in the system. These theories essentially treat the solid 
phase as an inhomogeneous liquid, in which the density 
profile consists of an array of density peaks, each corre- 
sponding to a localized particle, in contrast to the liquid 
phase which has a uniform density distribution. The PFC 
is a theory for an order parameter profile </>(x), which in a 
similar manner takes a constant value in the liquid phase 
and forms an array of peaks in the solid phase. 

In this paper we consider a simple liquid that has been 
rapidly quenched to a temperature well below Tf and de- 
velop a theory for how the solidification front propagates 
into the unstable liquid. We base our analysis on the 
DDFT and PFC models. DDFT predicts that the time 
evolution of the one-body density /5(x, t) of a system of 
particles is governed by the following equation [9HT2]: 



dp(x, t) 
dt 



= rv • 



p(x,i)V 



SF[p] 



(1) 



where T is a (constant) mobility coefficient and F[p(x, t)] 
is the equilibrium fluid Helmholtz free energy functional: 



F[p(x,<)] = /T 1 J dxp(x,t)[ln(p(x,i)A 3 )-l] 
+ F ex [p(x,t)] + J dxV ext (x,t)p(x,t). 



(2) 



The first term is the ideal gas free energy; A is the ther- 
mal de Broglie wavelength and /3 = 1/ksT is the in- 
verse temperature. The second term F ex is the excess 
contribution and the final term is the contribution from 
the external potential V cyi t (x, t) . The DDFT may be de- 
rived from the Smoluchowski (Fokker-Planck) equation 
for a system of interacting Brownian particles with over- 
damped stochastic equations of motion, by assuming that 
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the two-body correlations in the non-equilibrium fluid are 
the same as those in an equilibrium fluid with the same 
one-body density profile |TT] . Moreover, for dense atomic 
or molecular fluids, in which the equations of motion for 
the particles are, of course, Newton's equations of mo- 
tion, one can argue ES] EI] that Eqs. Q and ^ still 
provide a reasonable approximation for the dynamics of 
the system, particularly when it is not too far from equi- 
librium. 

This paper is laid out as follows: In Sec. [IT] we con- 
sider the stability of a uniform liquid with number den- 
sity p(x) = po and obtain the dispersion relation for the 
growth/decay of small amplitude harmonic density per- 
turbations. We then approximate this relation and ob- 
tain a simple expression which coincides with the disper- 
sion relation that one obtains from considering the PFC 
theory. In Sec. |III| we employ the marginal stability 
hypothesis to compute from this dispersion relation the 
speed of a solidification front advancing into a linearly 
unstable uniform liquid EI] ■ We also make an expansion 
in a certain small parameter related to undercooling, in 
order to obtain an analytical expression for the speed c of 
the solidification front. We find that the wavelength A of 
the density modulations which develop in the system as 
the solidification front advances, is not necessarily equal 
to the lattice spacing A c of the equilibrium crystal that 
the system seeks to form. This is because the length A is a 
dynamically selected (non-equilibrium) quantity. When 
the liquid is only weakly supercooled into the linearly 
unstable region, then A « A c and one should expect a 
regular crystal to form easily. However, when the system 
is deeply supercooled, then A ^ A c and one should ex- 
pect the formation of a regular crystal to be frustrated 
and the structure that is initially formed behind the ad- 
vancing solidification front to be somewhat disordered 
(amorphous). In Sec. IV we confirm this conclusion, i.e., 
that a deep quench leads initially to the formation of 
solids with greater disorder, using numerical simulations 
of the PFC model in two spatial dimensions. We also 
show that the transverse filamcntation of the stripe pat- 
tern nucleated by the advancing front is a consequence 
of the random initial conditions we employ. In Sec. |Vj 
we draw our conclusions and discuss the applicability of 
our PFC results to understanding real materials. 



II. DISPERSION RELATION 



We consider a bulk fluid, where the external potential 
14 x t (x, t) — in Eq. ([2| and consider small density fluc- 
tuations p(x, t) = p(x, t)—po about the bulk fluid density 
Po- We have in mind that we are considering a homoge- 
neous fluid which has been rapidly quenched to the region 
of the phase diagram where the crystal is the equilibrium 
phase. In the following derivation of the dispersion re- 
lation for the growth/decay of harmonic density fluctua- 
tions p(x, i), we initially follow Ref. [IT]. From Eqs. 



and ^ we obtain: 

10p(x,t) = v2 ^, ) _ /9oV 2 c (l) (X;i) 

-V.[p(x,t)Vc«(x,t)], 



D dt 



(3) 



where the diffusion coefficient D — T//3 and c^^x, t) — 
—f36F cx /5p is the one-body direct correlation function 
[5] IS]. We linearise Eq. ^ in p by Taylor expanding 
about the bulk fluid value, giving 



«(x) = c rc 



(oo) + 



, 5c«(x) 



<5p(x') 



p(x',t) + <D{~p 2 ) 



I'd 



( 4 ) 

where c^(oo) = c^'[po] = — /3/i ox and p ox is the excess 
chemical potential. Note also that 



fcW(x) 
Jp(x') 



= -P 
= c< 2 > 



S 2 F cx [p] 



(x,x') 
|x - x'|;p ) 



(5) 



for a homogeneous fluid of spherically symmetric parti- 
cles. For an equilibrium system c( 2 )(|x — x'|;po) is the 
Ornstein-Zernike direct pair correlation function of the 
fluid with density p . Substituting Eq. Q into Eq. (pjj), 
we obtain 



1 dp(x,t) 
D dt 

-PoV 2 



V 2 p(x,t) 



J dx' C ( 2 >(|x-x'|;p )p(x',t) 



0(p 2 ). (6) 



We now assume that the density fluctuation is of the form 
p(x, t) = eexp(Lut + ik.x), where e is a small amplitude 
and the dispersion relation w(fc), where k = |k|, is yet to 
be determined. From Eq. ^ we obtain: 



D 



p(x,t) = ~k 2 p{ x ,t) + p k 2 c(k)p(x,t) + O(p 2 ), (7) 



where c(k) = J dxexp(— ik.x)c^ 2 ^(x; po) is the Fourier 
transform of the pair direct correlation function. Note 
that for an equilibrium fluid, at a state point outside the 
spinodal, S(k) = (1 — p c(/c)) _1 is the static structure 
factor. Linearising Eq. ^ we obtain the dispersion rela- 
tion: 



}{k) = -Dk 2 [l - p c(k)}. 



(8) 



It is clear that small density fluctuations only grow in 
amplitude if for some wave numbers k we have u(k) > 
0. Crystallisation occurs when the system is unstable 
against periodic density modulations, which occurs when 
w(fc) > for a band of wave numbers about k as q, where 
q =/= 0. The dispersion relation w(fc) for an unstable sys- 
tem is of the form sketched using the solid line in Fig. 
[JJ Note that crystals may form before the system be- 
comes linearly unstable; however, in this case the crystal 
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FIG. 1: Sketch of the approximate dispersion relation u(k) 
in Eq. fl3b. When the dispersion relation takes the form 
labelled 'stable' the uniform liquid is linearly stable (A > 
0, dashed line). In the case labelled 'unstable' the uniform 
fluid is linearly unstable and density modulations with wave 
number k w q grow in amplitude, leading to the formation of 
the solid phase (A < 0, solid line). 



must be nucleated. Furthermore, if the fluid state falls 
within the solid-liquid coexistence region, then the crys- 
tal front will not advance indefinitely: it will grow until 
it has removed sufficient material from the surrounding 
fluid to produce phase coexistence between the liquid and 
the crystal. In this case, the crystal forms a 'localised 
state'. PFC results for this situation may be seen in 
Refs. 0H1 122]. 

In the following we assume that the speed with which 
the crystallisation front advances into the unstable liquid 
corresponds to the marginal stability criterion [25H27] . 
Specifically, we suppose that the unstable liquid state 
is characterized by a dispersion relation to = uj(k). In 
the frame in which the front is stationary, the dispersion 
relation becomes u) = ick + u(k) = Q(fc), where c is the 
speed of the front. In this frame the following relations 
hold 



= 



dk 

»(n) = o 



(9) 
(10) 



corresponding to the presence of a double root of u) = 
£l(k) in the complex fc plane together with the require- 
ment that the perturbation neither grows nor decays. 
Since 3(0) ^ in general, the wavetrain left behind 
by the moving front has a well-defined frequency in the 
frame of the front. 

The above conditions are equivalent to three conditions 
which are to be solved for the speed c of the crystalli- 
sation front together with the associated complex wave 
number k = k r + iki. The resulting density profile at or 
before the front has the form p(x, t) = /5f ron t(£, i), where 
£ = x — ct represents the position relative to the moving 
front and pi TO nt(^,t) ~ exp(— fcj£) sin(fc r £ + 3(0)i). Thus 



k r is the wave number of the growing perturbation, i.e., 
the wave number before the front, while ki represents 
the spatial decay (growth) of the perturbation in the for- 
ward (backward) direction. In contrast, the pattern left 
behind by the front is a fully nonlinear periodic state with 
wave number k* , say. In the absence of phase slips such a 
state takes the form p(k*^+^s(Cl)t), i.e., a wave that trav- 
els backwards relative to the front with frequency 3(0); 
with no phase slips this frequency is identical to the fre- 
quency ahead of the front and so can be computed from 
the marginal stability calculation. In view of the gradient 
structure of Eq. ([I]) this solution must be stationary in 
the laboratory frame. Thus p(k*x — k* ct + 3(Q)t) must 
be independent of the time t, implying that [26] 



k* = -9f(f2) = k r + 
c c 



(11) 



This equation expresses the conservation of nodes. Note 
that k* differs in general from the marginal stability wave 
number k r . 

To obtain the crystallisation front speed c, one assumes 
an approximation for F ex in Eq. ^ to obtain an expres- 
sion for c(fc) |42j . and hence the approximate dispersion 
relation w(fc). With this input Eqs. ^ and (10) may 
be solved (numerically) for c, k r and fcj and the wave 



number of the deposited pattern evaluated using (111. 
Under certain conditions an approximate solution to this 
problem may be obtained analytically, as shown next. 



III. SOLIDIFICATION FRONT SPEED 
A. Approximate dispersion relation 

To compute the front speed we first derive an approx- 
imation to the dispersion relation by expanding c(k) in 
powers of k. In order to capture the peak at k s=s q in the 
dispersion relation, one must retain at least terms up to 
0(fc 4 ) in c(k). Thus we write 



c(k) ~ Co + c-ik 2 + Cik 4 



(12) 



and suppose that C4 < 0. This approximation corre- 
sponds to making a gradient expansion of the free energy 
F C y\p] and retaining only terms up to and including the 

terms [V 2 p(x,i)] 2 43. Substituting Eq. ^ into 

Eq. {§), we obtain 



;(fc) = -afc 2 [A+(q 2 -/c 2 ) 2 ], 



(13) 



where a — —poc^D, q 2 — — C2/2C4 and A = (pqCq — 
1) / P0C4 — (C2/2C4) 2 . The uniform fluid thus becomes lin- 
early unstable for A < 0; i.e., the stable dispersion curve 
in Fig. [T] corresponds to a case when A > and the un- 
stable curve is for A < 0. Thus the magnitude of the 
parameter A indicates how deep one has quenched into 
the region of the phase diagram where the uniform liquid 
is linearly unstable. 
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Note that the dispersion relation in Eq. ( 13 ) is exactly 



that which one obtains when considering the PFC model 
for the order parameter ^(x, t) = [p(x, t) — po]/pi, where 
pi is a constant. The PFC model may be derived from 
the DDFT by assuming a gradient expansion in F ex and 
expanding the free energy in powers of 4> an d then lin- 
earising certain terms [T31 [151 HE] > obtaining 



dt a 6(f>(x,t)' 



(14) 



Here a is a mobility coefficient and the free energy func- 
tional 



F[4>] = J dx 



(15) 



Details of this derivation are contained in Appendix A. 
For the PFC model, we find that the uniform state 
</>(x, t) = 4>q (corresponding to the liquid) is linearly 
unstable when the undercooling parameter r < — 3(^q. 
Thus, in this model we have A = r + 3<j>Q and A < 
represents the undercooled liquid state. 

For the one-dimensional PFC model the marginal sta- 
bility analysis described above was performed in Ref. |21j . 
In the remainder of this paper we extend the predictions 
of this approach both analytically and numerically, and 
compare them with results from numerical simulations in 
one and two dimensions. 



for the wave number behind the front. 

These results show that when A is small, the front 
propagation speed c oc V— A. One also sees that fej cx 
y/—A, a cx y— A while k* — k r cx |A| and so increases as 
—A increases. The above results are accurate when |A| 
is small, but are not reliable when the system is deeply 
quenched, i.e., when | A| is not small. In particular, when 
this is the case, it is important to distinguish between 
the wave number k r predicted by the marginal stability 
condition and the wave number k* left behind by the 
moving front. In fact, one can obtain an expression for 
the crystallisation front speed c that is more accurate for 
a larger range of values of A as follows. We start by 
linearising the real part of Eq. ^ in ki to obtain 



qA 



o(A + 4g 4 ) 



(20) 



We next expand the imaginary part of Eq. ([9]) to second 
order in fcj and use Eq. ( 20 ) to obtain 



2aAq(16q s - 28g 4 A + A 2 
c(16<z 8 + 8g 4 A + A 2 ) 



(21) 



Together these results determine an approximation for 
k r = q + aki. Finally, we expand Eq. ( 10 1 to second order 
in fcj. Using Eqs. (201 and (21) leads to the following 



expression for the crystallisation front speed 



B. The front speed 

We now assume A < and calculate the speed with 
which the solidification front propagates into the unstable 
liquid. Taking the approximate dispersion relation in Eq. 
(13) together with Eqs. ^ and (10), hereafter the full 



theory, we obtain three equations for the three unknowns 
c, k r and ki |21j . Two of the resulting equations are 
quintics in k r and ki and the third has a term in k®. 
These simultaneous equations may be solved numerically. 
Results for the front speed c obtained from doing this are 
displayed in Fig. [2] (a) as a solid line. However, one can 
proceed further analytically by noting that when A is 
small ki is also small. We also make the ansatz that k r w 
q + aki, where the constant a is a variable to be solved 
for. We now proceed by expanding the three equations 
we obtain from Eqs. (|9| and ( 10 ) in powers of fej 



Oi 



can linearise all three equations in ki and then solve for 
c, ki and a to obtain the following: 



c = aq^/SAq 4 + 2A 2 
9v /-8Ag 4 + 2A 2 



hi — 



2(4g 4 + A) 
2A 

V-SAg 4 + 2/P 



(16) 
(17) 

(18) 



In addition, expanding Eq. (11) yields the prediction 
2a 

k* = k r qkAA + 4q 3 ah - 6q 2 fc 2 ) (19) 

c 



_ 4aq 3 (l6q 8 - 28<? 4 A + A 2 )^-A(A 2 - 64g 4 A + I6q 8 ) 
(A 2 - 64g 4 A + 16(f) (4g 4 + A) 

(22) 

C. Comparison with numerical simulations 

In Fig. [2ja) we display the result from Eq. |22| ) as 
the dashed line, together with the result from the full 
numerical solution (solid line) obtained from Eqs. j9j, 
( flO| and Q. We see that for small values of |A| < 0.1 
the expression in Eq. ( [22] ) for c is accurate. However, 
for larger values of |A[, it becomes less reliable. This 
approach also captures fairly well the behaviour of ki, as 
can be seen in Fig. [2jc) . However, as can be seen in (b) , 
it does not describe very well the behaviour of k r as a 
function of A. 

In Fig.[2|a) and (b) we also display results for the front 
speed c and the wave number k r obtained numerically 
by solving the PFC equations Eqs. (14) and (15) on a 



1-dimensional grid. We set the system size to be > 1000, 
which is sufficient for a stationary advancing front to de- 
velop [U]. We compare results obtained for various val- 
ues of the spatial grid spacing dx. We present results 
for dx = 0.2, 0.5, 1 and 7r/3 (in the literature there are 
some groups that use this particular value) . We find that 
for the larger values of the lattice spacing dx the front 
speed c is markedly slower than for smaller values of the 
lattice spacing, which are in good agreement with the 
exact speed obtained by solving Eqs. (13), ^ and (10) 
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FIG. 2: (Colour online) (a) The crystallization front speed c as a function of |A|, for q = 1, a = 1. The solid line is the 
result from solving the full theory [Eqs. Q, ( |10[ ) and 113k], and the dashed line is the analytical approximation in Eq. fl22p. 
Note that in the PFC model A = r + 3(f>Q. The symbols correspond to numerical results obtained for the front speed when the 
PFC Eqs. ( 14 1 and ( 15 I are discretized with various values for the spatial grid spacing dx, as indicated in the figure. In (b) we 
display the corresponding results for k r and in (c) for ki. In the inset of (c) we display the order parameter profile plotted as 
log \ 4>(x) — 0o | versus x, for the case when (fio = —0.4 and r = —0.9, i.e., A = —0.42. We did not extract the value of ki from 
the numerical results. However, as can be seen from the inset to (c), the agreement between the slope of the dashed line, which 
has gradient — fe, as obtained from the full theory, and the envelope of the order parameter profile, is good. Panel (d) shows the 
corresponding results for the wave number k* of the nonlinear state deposited by the front (symbols) for comparison with the 
predicted wave number k* (solid line). The inset shows a plot of both k* and k r which confirms that the wave number behind 
the front differs from k r , the wave number amplified by the front (dashed line). Both k* and k r differ substantially from q = 1. 



numerically and displayed as the solid line in Fig. [2ja). 
Finally, in Fig. [2jci) we display the corresponding results 
for k* and compare these with the theoretical predictions 
for k* (solid line) and k r (dashed line). The theoretical 
predictions for c, k r and k* (solid lines in Figs. [2]ja,b,d)) 
are in excellent agreement with the numerical results ob- 
tained with grid spacing dx — 0.2. Figure [2] also shows 
that results obtained with dx > 0.5 are substantially in 
error. This is because the discretisation of the system 
effectively adds a friction term proportional to the mag- 
nitude of dx to the dynamical equations, which slows 
down the advancing front - i.e., the numerical grid can 
'pin' the advancing front. Evidently, this pinning effect 
is also reflected in the corresponding values of k r and k* . 

We did not extract the value of ki from the numerical 
results. However, in the inset of Fig. [2jc) we display the 
order parameter profile plotted as logics) — </>o| versus 
x for the case when O == —0.4 and r = —0.9, calculated 
numerically using the grid spacing dx = 0.2. From the 



analysis of the advancing front profile one expects that 
4>{x,t) — (f>o ~ exp(— ki(x — ct))sm(k r x), so that when 
the order parameter profile is plotted in this manner, 
the envelope function exp(— kix) of the advancing front 
profile becomes a straight line with gradient — ki. The 
dashed line in the inset of Fig.[2|c) is a straight line with 
gradient — ki computed from Eqs. ( fl3| ), ^ and (10). It is 
clear that the gradient of the envelope of the numerically 
obtained order parameter profile is very close to that of 
the dashed line. Thus, we conclude that the analysis 
based on Eqs. (13 1, Q and (10) leads to a prediction 
for the solidification front speed c which is precisely that 
which one obtains from solving the PFC equations ( 14 ) 
and ([15]). 

In Fig. [3] we show an example of a front propagating 
towards the right when A = —0.42. The front region 
is clearly visible on the semi-logarithmic plot shown in 
panel (b); panels (c) and (d) show enlargements corre- 
sponding to the region behind the front and the front 
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FIG. 3: (Colour online) (a) A front advancing to the right at one instant of time when r = —0.9, (j>o = —0.4, q — 1 and a = 1, 
computed with dx = 0.2. (b) The solution in panel (a) on a semi-logarithmic plot, (c) Enlargement of the region behind the 
front in panel (a), (d) Enlargement of the front region in panel (b). 



region itself. From these figures one determines that the 
wave number in the front region is k r « 1.189, while the 
wave number behind the front is k* s» 1.123. These mea- 
surements agree very well with the exact marginal sta- 
bility result, k r ss 1.187, and the prediction in Eq. (Ill, 
k* « 1.129. 

It is of interest to note that the dynamically selected 
wave number k*, which determines the wavelength A = 
2-w/k* of the density modulations left behind the advanc- 
ing front, can differ significantly from the equilibrium 
wavelength A c w 27r/<7 of the fully formed crystal. This 
means that for large negative values of A, which cor- 
responds to a deep quench (i.e., the unstable liquid is 
strongly supercooled), the system must perform signifi- 
cant rearrangements after the initial solidification front 
has passed, in order to obtain density modulations with 
wavelength 27r/g, corresponding to an ordered crys- 
tal of minimal energy. However, one should expect that 
these later rearrangements (ageing) are frustrated by the 
fact that the system has already chosen a different and 
dynamically selected length scale. As the system ages 
some of these defects anneal, reducing the disorder in 
the solid and bringing it closer to equilibrium. We 
thus believe that the difference between the dynamic and 
equilibrium crystalline wavelengths may be an important 
factor in understanding why some rapidly quenched liq- 



uids and soft matter systems exhibit disorder rather than 
forming a regular crystalline material. We illustrate and 
demonstrate this result further in the next section. 



IV. PFC RESULTS 

In this section we confirm and illustrate the results and 
conclusions of the analysis given in the previous section, 
using results obtained from direct numerical simulations 
in two spatial dimensions for the simple PFC model given 
in Eqs. (|14[) and (151. This model system is now well un- 



derstood and much is known about its thermodynamics 
and phase behaviour, and the structures that are formed 
|H [T3HH]. We display in Fig. [I] the order parameter 
profiles for a solidification front advancing from left to 
right into the unstable uniform liquid phase, for a system 
with q = 1, (f)Q = —0.43 and r — —0.9, corresponding to 
A = —0.35. The profiles in Fig. [4] are calculated by tak- 
ing an initially uniform system with 4>{x) = 4>q, of size 
2000 x 50 with periodic boundary conditions and grid 
spacing dx — dy = 0.5. The solidification is initiated 
by adding small amplitude random noise to the profile 
along the line x = at time t = 0. The order param- 
eter profiles displayed in Fig. [4] correspond to the times 
t* = at/q 2 = 140, 200, 260 and 340. We see that the 




50 




300 



350 



400 



450 



500 



550 



600 



650 



700 



FIG. 4: (Colour online) Order parameter profile s fo r a crystallisation (solidification) front advancing into the unstable 
uniform phase, obtained from the PFC model (Eqs. (14 1 and (15}) in two spatial dimensions when cj>o = —0.43 and r = —0.9, 
corresponding to A = —0.35. The plots correspond to the times t = at/q 2 — 140, 200, 260 and 340, going from top to bottom. 
The solidification front was initiated at t = at x = and propagates towards the right. Note the rearrangements that occur 
at points well behind the moving front. 



front advances by first forming stripc-likc density mod- 
ulations in the direction of travel, as predicted by the 
analysis in Sec. |III| above. However, the stripes are typ- 
ically broken into transverse domains or 'filaments' (see 
Sec. IV A I, leading to a two-dimensional structure that 



subsequently breaks up into density peaks resembling a 
solid. Figure [4] shows the order parameter profile corre- 
sponding to the time t* — 260, and reveals that there is 
a significant amount of disorder in the arrangements of 
the density peaks shortly after the solidification front has 
passed. Then, over time, the system rearranges (ageing) 
leading to the more regular ordering seen in the order 
parameter profile for the time t* = 340 (see Sec. IV B I. 



The results of the PFC model depend on both the cho- 
sen value of A < 0, the undercooling, and of <fio, the 
background homogeneous state into which the solidifica- 
tion front propagates. To explore the parameter space, 
we solved the PFC model of size 400 x 400 with periodic 
boundary conditions in the y direction and dx — dy — 
0.5, initializing the solidification front by adding small 
amplitude random noise to the order parameter profile 



along the line x = at the time t = [35] • We use the 
same realization of the initial condition throughout. Ac- 
cording to the theory presented in the previous section, 
the front speed c and wavenumbers k r and k t are deter- 
mined by the value of A only. Figure [5] shows the results 
for A = —0.1 and several different values of r (equiva- 
lently of 4>o, since 4>o = \/(A — r)/3), all at the same time 
t* = 152 after the front was initiated at x = at t = 0. 
Detailed analysis shows that the front speed and length 
scale right in the front region are indeed independent of 
the background value of (f>o- On the other hand, the ex- 
tent of the region of the stripe-like state is dramatically 
reduced as |<^>o| increases. Since mature stripes of wave- 
length 27r/fc* are created at the rate 3(f2) and destroyed 
at the rate Whex at which the instability to hexagons man- 
ifests itself, it follows that the width £ of the stripe region 
scales as £ ~ 2%^s(XY)/k*ui} lex — 27rc/u;i iex using Eq. (Ill, 

cf. E 

like 



J1I30]. As shown in Appendix B, this quantity scales 
6o | 1 with a coefficient of proportionality that is in- 
dependent of A when |A| <C 1. Our numerical results are 
consistent with this prediction, although it is somewhat 




FIG. 5: (Colour online) Order parameter profiles of a crystallisation (solidification) front advancing into the unstable uniform 
phase, obtained from the PFC model (Eqs. ( |14[ ) and ( |15[ ) ) in two spatial dimensions for A = —0.1 and several different values 
of r, all at time t* = at/q 2 = 152. The solidification front was initiated at x = at t = and propagates towards the right. 
The top left panel is for r = —0.2 and <j>o = —0.183, top right for r — —0.5 and (j>o — —0.365, bottom left for r = —0.9 and 
<j>o — —0.516 and bottom right for r = —1.3 and </>o = —0.632. The displayed region is part of a larger system of size 400 x 400 
and calculated with a grid spacing dx — dy — 0.5. Note the different colour table scales in each panel. 



difficult to determine precisely the width £ from the data. 
In fact, simulations starting from random initial condi- 
tions show that for low |<^o| the instability of the stripe 
state generates structures that are more rhomboid than 
hexagonal. With increasing |</>o| the structures become 
more hexagonal but the fraction of vacancies within the 
structure goes up. This is a consequence of the fact that 



the curve A(0o) = —0.1 in the (</>o, r ) plane moves as c/)q 
increases and eventually crosses into the coexistence re- 
gion between the hexagonal crystal and the homogeneous 
or liquid state [22] . 
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A. The transverse length scale 

Figures [4] and [5] reveal the presence of unambiguous 
filamentation of the stripe pattern created by the pas- 
sage of the front. To understand the origin of this fil- 
amentation we show in Fig. |6][a) the quantity </>(x) at 
time t* = at/q 2 = 64 when A = -0.2 (0 O = -0.483 
and r = —0.9) while Fig. [6jb) shows the same solu- 
tion but in terms of the quantity In |</>(x) — 4>q\. The 
latter representation not only rectifies the solution, but 
also amplifies it strongly in regions where </>(x) 4>q. 
The figure reveals that the filamentation is present al- 
ready at the front where the amplitude of the stripes is 
still minute, of order e -16 . Careful study of the origin 
of this filamentation shows that it is a consequence of 
the perturbation used to initialize the simulation. The 
ridges that break up the stripe pattern correspond to 
zero-crossings in <f>(x = 0, y, t = 0) — 4>q, here a particular 
realization of a uniformly distributed random variable 
on the interval [(f> — 0.1, O + 0.1]. The regions where 
<f)(x = 0, y, t = 0) — 4>q w travel more slowly than re- 
gions where (j)(x = 0, y, t = 0) — <j)o 7^ 0, and the latter are 
broad enough to trigger the formation of stripe segments. 
Thus the filaments are an imprint of the initial condition, 
and the advancing front acts as a noise amplifier. Sim- 
ulations initialized from a small amplitude perturbation 
with a single wavenumber k± preserve this wavenumber 
into the nonlinear regime and the resulting filamentation 
is periodic with wavelength Aj^ = 2n/k±. 

The stripes created by the advancing front are unsta- 
ble to oblique disturbances that favour the formation of 
hexagonal structures and this instability becomes visible 
once (j>— (f>Q = O(4>o) (Fig.[6]ja,b)). The growth rate of this 
instability is proportional to \(f>o\ (see Appendix B) and 
consequently we expect the width t of the stripe interval 
ahead of the hexagonal pattern to decrease with increas- 
ing | 4>q |, all other parameters remaining fixed (cf. Fig.[|). 
Our simulations reveal, however, that the filamentation 
imprinted by the initial conditions also has a strong effect 
on the ability of the system to form hexagons. If the char- 
acteristic transverse scale Aj_ is far from 2Ay/v3, where 
An = 2ir/k*, we find that the formation of hexagons is 
delayed until such time as the required wavenumber is 
generated by nonlinear interactions. Thus the initial con- 
dition strongly influences, through the above process, the 
time required to form the crystalline state. Moreover, 
since the selected wavelength Aii is likewise non-optimal, 
both factors contribute to frustration and disorder in the 
solidification process for deep quenches. 



B. Structure and correlations over time: ageing 

In order to quantify the degree of order in the system 
and compare results from shortly after the solidification 
front has passed with those at a later time, we computed 
the bond angle distribution p(8) and radial distribution 
function g(r) as a function of time after the solidification 
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FIG. 6: (Colour online) Order parameter profile at time 
i* = at/q 2 — 64 for a crystallisation (solidification) front 
advancing into the unstable uniform phase, obtained from 
the PFC model (Eqs. ( |14[ | and ( |15[ )) in two spatial dimensions 
when (j>o = —0.483 and r = —0.9, corresponding to A = —0.2. 
The solidification front was initiated at x = at t = and 
propagates towards the right. The upper figure shows the 
order parameter profile <^(x), while the lower figure shows the 
same profile, but instead plotting the quantity In j</>(x) — <^o|- 
Plotting this quantity reveals the fine structure in the profile 
ahead of the front - note the scale: the smallest amplitude 
structures that are displayed have an amplitude ~ e -16 ). The 
displayed region is part of a larger system of size 400 x 400 
and calculated with a grid spacing dx = dy — 0.5. 
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FIG. 7: (Colour online) The bond angle distribution p(9) (top row) and radial distribution function g(r) (bottom row) at 
various times t* , after the solidification front was initiated. The undercooling parameter is r = —0.9, and the value of <f>o 
decreases from left to right (as indicated in the figures) resulting in (from left to right) A = —0.42, —0.29 and —0.15. After the 
initial crystallisation front passes by, the system undergoes 'ageing' as the particles are able to perform some rearrangements. 



front was initiated. These quantities are calculated from 
larger scale (grid size 400 x 400) simulations by first lo- 
cating all the maxima in the order parameter profile, i.e., 
the coordinates of all the density peaks (particles) after 
the crystallisation front has moved through the system. 
From these sets of particle coordinates, we calculate the 
radial distribution function g(r) in the usual way |31j . 
Since g(r) is a spatial two-point correlation function, it 
gives the probability of finding another particle at a dis- 
tance r away from any other given particle |32j . The bond 
angle distribution function is calculated by performing a 
Delauney triangulation on the system. The histogram of 
the values of the corner angles of this set of triangles (i.e., 
the nearest neighbour bond angles) is p(0). 

In Figs.CTa)-(c) we display the bond angle distribution 
p(6) for r — —0.9 as it varies over time, for (a) 0o = —0.4, 
(b) 4>o = —0.45 and (c) <fio = —0.5. These three values 
of 0o correspond, respectively, to A = —0.42, —0.29 and 
—0.15. These bond angle distributions are centred on 
the value 60°, due to the dominant hexagonal ordering 
in the system, and we see no peaks at 45° and 90°, which 
would indicate square ordering [22]. We see in (a) and 
(b), corresponding to larger values of |A| (i.e., the deeper 
quenches), that at the time t* — 200 the distribution p{6) 
is much broader than for later times, indicating that at 
this early time, shortly after the solidification front has 
passed through the system, there is much more disorder 
in the system than at the later times. Over time, the sys- 
tem rearranges to form a much more ordered solid, with 
p(8) being much more sharply distributed around 60°. In 
contrast, for the shallow quench case with small |A| dis- 



played in Fig.JTJc), we see that p(9) is sharply distributed 
around 60° even for short times after the solidification 
front has moved through the system and that it does not 
change much as time goes by, indicating there is very lit- 
tle ageing in the system. These findings can also be seen 
by inspecting the radial distribution functions g{r) dis- 
played in Figs. |7jd)-(f). For the shallow quench case in 
(f) we see that g(r) does not change much over time. In 
contrast, for the deeper quench cases in (d) and (e) we see 
that at t* = 200 the decay g(r) — 5- 1 is much faster than 
at later times. The fact that the amplitude of the oscil- 
lations in g(r) is much smaller at earlier times indicates 
that there is much less long range (crystalline) ordering 
in the system. As time proceeds, the amplitude of the 
oscillations in the tail of g(r) grows, indicating that the 
system is rearranging to form a much more ordered sys- 
tem with the particle locations being well correlated over 
larger distances. The larger amount of disorder shortly 
after a deep quench is a consequence of the mismatch 
between the wavelength selected dynamically by the ad- 
vancing solidification front and the equilibrium lattice 
spacing of the crystalline solid. This mismatch increases 
with increasing |A|. The initial appearance of density 
modulation with the 'wrong' wavelength creates disorder 
and frustration in the system, a picture corroborated by 
the results in Sec. Mil 
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V. CONCLUDING REMARKS 

In this paper we have studied the propagation of a 
solidification front into a supercooled liquid, i.e., into 
a linearly unstable state. We employed dynamical den- 
sity functional theory (DDFT) to derive an approximate 
dispersion relation for small perturbations of the spa- 
tially uniform liquid state and noted that this disper- 
sion relation is identical in form to that derived from 
the phase field crystal (PFC) model of crystal growth. 
In both approaches the solid phase is represented as 
a spatially structured state with local maxima in the 
density profile p(x) or equivalently the order parameter 
</>(x) representing the time-averaged location of individ- 
ual atoms/particles. The present approach is thus able 
to bridge purely continuum or macroscopic solidification 
theory |33) with atomistic approaches such as molecu- 
lar dynamics. Despite this fundamental difference, the 
DDFT and PFC models that result can still be formu- 
lated in terms of partial differential equations. These 
may be nonlocal as in DDFT or local as in PFC. 

Knowledge of the dispersion relation suffices for the 
computation of the speed of the solidification front when 
this speed is selected by linear processes, i.e., in situa- 
tions where the growth of the perturbations behind the 
front compensates for the propagation of the front, re- 
sulting in a steadily advancing front of constant shape. 
However, in some problems the speed of the front may 
instead be determined by nonlinear processes pH]. For 
this reason it is essential to compare the prediction ob- 
tained from the linear marginal stability criterion em- 
ployed here in the form of Eqs. ([£])-( 10 ) with numeri- 
cal simulations. Such simulations yield in addition im- 
portant information about processes occurring on longer 
time scales than the propagation time. Our results can 
be summarized as follows. For small undercooling, as 
measured by the parameter |A|, the advancing front se- 
lects wavelengths close to the equilibrium wavelength A c 
of the crystalline solid, resulting in steady transforma- 
tion of the liquid state into solid. The front speed is 
c ~ \/~A. For large undercooling (i.e., supercooling) 
the front speed is faster and follows the approximate re- 
lation c ~ — A. In this regime the wavelength selected 
by the advancing front differs substantially from A c re- 
sulting in a nonequilibrium structure that subsequently 
evolves on a longer time scale, first via an instability to 
a hexagonal structure and subsequently via slow defect 
migration and annihilation. This 'ageing' process con- 
sists of rearrangements as the system seeks to anneal out 
the defects and differently orientated domain structures 
which frustrate the formation of a regular crystal with 
wavelength A c . 

We have also found that the initial perturbation im- 
printed on the advancing front may have a significant 
effect on the manifestation of the instability of the stripe 
state with respect to hexagonal perturbations. Since this 
transverse scale will also differ from the optimal scale 
27r/ fc* its presence provides an additional source of frus- 



tration following the passage of the front. Although these 
results were obtained using the PFC model, analogous 
two-dimensional calculations based on a DDFT model 
yield very similar results (not shown). 

An important issue that we must mention concerns the 
extent to which insights from the PFC model can be ap- 
plied to solidification in real materials. The transition 
in the PFC model from uniform to modulated phase is 
weakly first order, stemming from a truncated gradient 
expansion approximation to obtain the PFC free energy 
functional in Eq. ( 15 1 (see also Eq. ( A6 1 in Appendix A). 
The fact that for some values of cf>o the PFC model ex- 
hibits a stripe phase that is not seen in real atomic fluids 
is an indication that the truncated gradient expansion 
approximation has failed for these 4>q values [22] . Thus, 
great caution should be taken in relating our results to so- 
lidification and glass formation in quenched liquids. For 
understanding how fronts propagate into a linearly un- 
stable fluid, the approach described above appears to be 
valid. However, owing to the very simple nature of the 
PFC model, we expect that its description of the struc- 
tures formed behind the front may be less reliable. The 
presence of a weakly first order transition in the PFC 
model makes it somewhat unrealistic as a model for ma- 
terials like liquid metals, but for soft matter (polymeric) 
systems we believe it is a good approximation. Much 
more work comparing the PFC to more sophisticated 
DDFT approaches, such as that presented in Ref. [IB] , 
is required in order to elucidate the extent to which the 
PFC can be used to model real materials. 

We mention, finally, that the DDFT presented above 
was derived for Brownian particles. Improvements in the 
theory required for application to atomistic fluids include 
the DDFT El El: 



d 2 p{^t) , Jp(x,t) 



1 



dt 2 



dt 



p(x,i)V 



SF[ P ] 
Sp(x,t) 



, (23) 



where m is the mass of the atoms and v is the collision 
frequency given by v « ksT/mD. Here D is the self- 
diffusion coefficient. The free energy F is given by 
Front propagation in one-dimensional models of this type 
is considered in EI]. Extensions of the present work to 
this class of models in two or more dimensions, together 
with a comparison with direct numerical simulations, will 
be presented elsewhere. 
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Appendix A: Derivation of the PFC model from 
DDFT 

The DDFT in Eqs. Q and ^ is a microscopic theory 
which describes the time evolution of the fluid one-body 
(number) density profile p(x, t) for a fluid of Brownian 
particles. In this section we start from the DDFT to 
derive the PFC model in its commonly used form. In 
our derivation we closely follow the arguments laid out 
in Ref. [13] . The excess contribution to the free energy 
_F CX in Eq. ^ is usually an unknown quantity. Here, 
we make an approximation for F ex , by making a Taylor 
series expansion in powers of p(x) = p(x) — po, where po 
is a reference density, giving [BJ: 



F ex [p(x)] = F ex [p ] + /dx p(x 



MM*)] 



\ J J dxdx! p(x)p(x' 



S 2 F ex [p(x)] 



2 . 



I'D 

(Al) 



The functional derivatives of the excess free energy which 
enter into Eq. (Al) are related to the n-body direct cor- 



relation functions in the following way [BJ: 
5 n F ex [p(x)} 



-k B Tc {n \ 



Xl,X 2 , 



<5p(xi)(5p(x 2 ) • • • <5p(x„) 

In particular, the first member of this series is the one- 
body direct correlation function, shown earlier in Eqs. ^ 
and Q. Note that the one body direct correlation func- 
tion evaluated in the bulk is equal to the excess chemical 
potential — fe^Tc' 1 ) (x) I = p ox . The second member of 



the series in Eq. ( A2 ) is the direct pair correlation func- 
tion: 



6 2 F c: 



6p(x)5p(x>) 



-fc B Tc (2) (x,x'). 



(A3) 



Substituting these expressions for the functional deriva- 
tives into Eq. ( Al ) and neglecting third and higher order 
terms we obtain: 

F ex [p(x)] fa F cx [p ] + Pox Jdx p(x) 

J J tfxdx' p(x) C ( 2 ) (x, x')p(x') (A4) 



The second term in this equation corresponds simply to 
a shift in the chemical potential and so this approxima- 
tion is commonly used without the second term explicitly 
written down 13, 35, 36 as originally done by Ramakr- 
ishnan and Yussouff [37] . To derive the PFC free energy, 
we make a gradient expansion of the two body direct cor- 
relation function and truncate at the fourth order term, 
giving [HIES]: 



,(2) 



(x, x') » -(3(C + C* 2 V 2 + C 4 V 4 )<5(x - x'), (A5) 



where in principle all the coefficients (7, are functions of 
p(x), although we assume here that the coefficients (7 2 
and C 4 are in fact constants. Inserting approximation 
((ABl into Eq. (lA4| gives: 



F ex [p(x)] « F ex [p ] + p c 



1 



Jdx p(x) 



+ ^ J dx p(x)(C + (7 2 V 2 + C 4 V 4 )p(x), (A6) 

which makes i 7 ' cx [p(x)] a local functional. Using this ex- 
pression for the excess free energy term we can now write 
the Helmholtz free energy for the system as: 



rfx 



F[P(X)] : 

where 
fo(p) = k B Tp(\n(p) 



/ (p(x)) + ^p(C* 2 V 2 + (7 4 V 4 )p 



(A7) 



1) + /cx[Po] + PcxP + ^C {p)p 



7 ;2 



(A8) 

Here the first term in f (p(x.)) comes from the ideal gas 
contribution (see Eq. ([2])) and we have assumed that the 
external potential V ex t — 0. We also make a further ap- 
proximation by making a Taylor expansion of the func- 
tion /o(p) around the reference density p , giving: 



Hp) 



fM + f'MP + ^-p 2 

, jf (P0) -3 ■ ftHpo) -* 

P At P ' 



3! 



4! 



(A9) 



We choose the reference density po so that the third 
derivative of the function /o(p) vanishes at p = po, 
i.e., /q (po) = 0. This gives the following: 



fo(p)~ fo(po) + fo(po)p+ 



fo(Po)~, , /o (4) (Po), 



-p 



4! 



p 4 . (A10) 



We now introduce a change of variables. We use the non- 
dimensional variable </> = p/pi, where p\ is a constant 



density, so Eqs. \K1\ and (|A10 1 become: 

1 



F 



(x)]= jdx 



where C 2 = C 2 /p 2 , C4 = C4/ ' p\ and 
fo(4>) ~a + l 



# 4 



(All) 
(A12) 



where a, b, c and d are constants. 

We now consider the dynamics of the model. We start 
with the DDFT equation ([!]). In the limit where p\(f> is 
small, the density preceding the gradient of the functional 
derivative becomes constant, i.e., p = po + pi4> ~ Po & n d 
Eq. ([I]) reduces to the following equation: 



<9p(x,t) 
dt 



TpoV 



2 5F[p(x,t)} 
Sp(x,t) 



(A13) 
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This is often referred to as "model B" dynamics in the 
classification of Hohenberg and Halperin [39]. Equiva- 
lently, we have the following equation for the time evolu- 
tion of the order parameter 0(x, t): 



d4>(x,t) = aV2 5F[<t>(x,t)} 



dt 



(x,t) 



(A14) 



where a = Tp / p\ is the mobility coefficient. Since the 
constant and linear terms in Eq. (A12) are irrelevant for 
the dynamics, we may drop the terms a + b(f> from the 
function /o(</>) in Eq. (A12|. The functional derivative of 



the free energy is then given by the expression: 



(A15) 



We may absorb the parameter d into the mobility co- 
efficient a. Also, we may choose p\ so that C±/d = 1. 
Writing C^/d = 2q 2 and c/d = r + q 4 , we arrive finally 
at the commonly used PFC free energy: 



F[0(x)] = /dx/(0(x)), 



(A16) 



where 



i 



2 T 4 



0(2<j 2 v 2 + v 4 ; 



(A17) 



Inserting these parameter values into the functional 
derivative of the free energy (Eq. (A15)), we obtain 
S -F = ( r + q 4 )<j) + <fi + 2g 2 V 2 + V 4 0. The PFC model 



is then given by the conserved dynamics in Eq. (A14|, 
where the free energy is given by Eqs. (|A16| and {Km. 



Appendix B: Instability of the stripe state 

In this Appendix we determine the timescale of the 
instability of the stripe state. This instability leads to the 
formation of the hexagonal structures shown in Figs. [5] 
and|6] 

We write the PFC model in the form 

4> t = aV 2 [A4> + (q 2 + V 2 ) 2 4> + 30 o 2 + 3 ], (Bl) 
where A = r + 30g and 

4> = <P~ 4>a 

= Ae ikx + Be ik (~ x+ ^ y ^ 2 + Ce ik (~ x ~^ v ^ 2 
+c.c. + h.o.t., (B2) 

where c.c. denotes the complex conjugate of the preced- 
ing terms and h.o.t. denotes higher order terms. Here 
A is the small but complex amplitude of the longitudi- 
nal mode while B and C are the corresponding ampli- 
tudes of two symmetry-related oblique modes. The state 



(A, B, C) = (A, 0, 0) thus corresponds to the stripe state 
while (A,B,C) = (A, A, A) corresponds to the hexagon 
state, with A > representing a hexagonal array of spots 
and A < representing a hexagonal array of holes or va- 
cancies. 

Weakly nonlinear theory now leads to the following 
equations for the amplitudes A, B,C: 



A t = -ak 2 [AA + 6<j) BC 
B t = -ak 2 [AB + 6(f> CA- 
C t = -ak 2 [AC + 60o AB - 



(B3) 
(B4) 
(B5) 



where the overbar denotes complex conjugation and A = 
A+(<7 2 — A: 2 ) 2 represents the bifurcation parameter shifted 
in proportion to the departure of the wavenumber k away 
from its optimal value k = q. By applying appropriate 
translations we may take A, B, C to be real. We also 
take B = C in order to focus on the instability of the 
stripe state with respect to hexagon-forming perturba- 
tions. The linear instability of an (j4oj0,0) state with 
respect to such perturbations is then described by the 
equation 



Bt 



-ak z \A- 



i>oA ]B, 



(B6) 



implying that the growth rate Whe 
bility is given by 



of the hexagon insta- 



Whox = -ak 2 [A + 6<p A ]. 



(B7) 



Here Aq is the amplitude of the stripe state. Within 



Eq. (B3) this amplitude is not determined: the growing 
stripe state (A < 0) does not saturate. However, the 
saturation amplitude of the stripe phase can be computed 
by setting B = C = and extending the above approach 
to cubic order while imposing the requirement that (0) = 
0, where (• • • ) denotes an average over the domain. We 
obtain A% = -4A(3-2^/g 4 ) -1 . For A«l these results 
(with A replaced by A) apply to stripes with k — k* since 
k* « q. 

Since A < for instability of the liquid phase, and 
likewise 0o < 0, the growth rate Whex is positive for 
all Aq > with k* ~ q, implying that the stripe state 
is always unstable with respect to the formation of the 
hexagon state with A = B = C > 0, i.e., a hexagonal 
array of spots. In the case 4>q > 3q 4 /2 the bifurcation 
to stripes is subcritical and the hexagon instability then 
competes with an amplitude instability. However, near 
threshold k* w q and the growth rate of the latter is 
therefore 0(|A|) while the growth rate of the hexagon 
instability is 0(|</>o| y|A~|) and so is larger. In either case, 
the longitudinal width I of the band of stripes ahead of 
the hexagonal state is predicted to scale, for small |A|, 
as I ~ (27rc/w hox ) +7i ~ 70 1 | _1 + 7i> where 70 is in- 
dependent of |A| but 71 cx fc" 1 does depend on |A|. For 
larger |A| the approximation in Eq. (22) is useful. 

The saturated hexagon state can be included selfcon- 
sistently in the above theory only when |^>o| <C 1, i.e., 
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when A 2 ~ \<fio\A ~ A [301 HI]- This is not the case in our simulations and we do not pursue this approach. 
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